function [pX,pY,OD] = ProcessImageDoubleGauss(hdfFilePath,autoROI,userROI,doFitData)

% [pX,pY,OD] = ProcessImage(hdfFilePath,autoROI,userROI,doFitData)
% Finds Gaussian-shaped signals on fluorescence or absorption images.
% Returns:
% pX - The parameters for a Gaussian fit to the data in the x (vertical)
% direction.
% pX(1) = Gaussian amplitude
% pX(2) = Peak position (in units of image pixels)
% pX(3) = Standard deviation (in units of image pixels)
% pX(4) = Offset
% pY - Same description as pX
% OD - The image (OD if it's an absorption image, raw image for a fluor.).

expp = ExpParams();
cons = Constants();

y = double(h5read(hdfFilePath,'/Image/image'));

dim = size(y);
abs = 0;

% Filter image to find ROI
N = 60;
sig = 30;

if size(dim,2) == 3 % Absorption
    num = y(:,:,1) - y(:,:,3);
    den = y(:,:,2) - y(:,:,3);
    den(den==0) = 1;
    ODlin = (num./den)';
    
    OD = -log(ODlin);
    OD(isinf(OD)) = 0;
    OD(isnan(OD)) = 0;
    OD = real(OD);
    OD = OD - mean(mean(OD));
elseif size(dim,2) == 2 % Fluorescence
    OD = y';
    OD = OD - mean(mean(OD));
else
    error('Can''t determine imaging method');
end

ODfilt = conv2(OD,gaussian2d(N,sig),'same');

yC = 1:expp.h_px_nb;
xC = 1:expp.v_px_nb;

if autoROI
    % Find ROI and extract first approx. signal from filtered image.
    mx = max(max(ODfilt));
    [I,J] = find(ODfilt == mx);
    wndw = 50;
    
    ROI(1) = I(1);
    ROI(2) = J(1);
    
    yIntFilt = sum(ODfilt(:,ROI(2)-wndw:ROI(2)+wndw),2);
    xIntFilt = sum(ODfilt(ROI(1)-wndw:ROI(1)+wndw,:),1);
    
    ampGuessY = max(yIntFilt) - min(yIntFilt);
    ampGuessX = max(xIntFilt) - min(xIntFilt);
    
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pYtmp = fminsearch(@(p)FitDoubleGauss(p,yC,yIntFilt),[ampGuessY ROI(1) 100 min(yIntFilt)],options);
    pXtmp = fminsearch(@(p)FitDoubleGauss(p,xC,xIntFilt'),[ampGuessX ROI(2) 100 min(xIntFilt)],options);
    
    % Use data from filtered image to extract data from raw image.
    sigmas = 3;
    xWndw = round(sigmas*pXtmp(3));
    yWndw = round(sigmas*pYtmp(3));
    xROI = ROI(2) - xWndw:ROI(2) + xWndw;
    yROI = ROI(1) - yWndw:ROI(1) + yWndw;
    
    % Remove parts of the ROI if it goes outside of the photo.
    if min(xROI) < 1
        xROI = xROI(xROI > 0);
    end
    if min(yROI) < 1
        yROI = yROI(yROI > 0);
    end
    
    if max(yROI) > expp.h_px_nb
        yROI = yROI(yROI < expp.h_px_nb + 1);
    end
    if max(xROI) > expp.v_px_nb
        xROI = xROI(xROI < expp.v_px_nb + 1);
    end
    
    yInt = sum(OD(yROI,xROI),2);
    xInt = sum(OD(yROI,xROI),1);
    
else
    % Use user-input ROI parameters to guess fitting parameters
    yROI = userROI(3):userROI(4);
    xROI = userROI(1):userROI(2);
    
    yInt = sum(OD(yROI,xROI),2);
    xInt = sum(OD(yROI,xROI),1);
    
    pYtmp(1) = max(yInt) - min(yInt);
    pXtmp(1) = max(xInt) - min(xInt);
    
    pYtmp(2) = [mean(userROI(3:4))];
    pXtmp(2) = [mean(userROI(1:2))];
end

if doFitData
    
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pY = fminsearch(@(p)FitGauss(p,yROI,yInt),[pYtmp(1:2) 150 min(yInt)],options);
    pX = fminsearch(@(p)FitGauss(p,xROI,xInt'),[pXtmp(1:2) 150 min(xInt)],options);
    
    % Algorithm:
    % - fit big Gaussian
    % - find indices of points above half
    % - fit Gaussian to this
    % - subtract Gaussian
    % - find position of mini-Gaussian
    % - fit double Gaussian to this, forcing mini-Gaussian to be at the
    % position of the trough in the subtracted signal.
    
    miniXint = xInt - gauss(pX,xROI);
    [miniGaussVal,miniGaussPosInd] = min(miniXint);
    mnGaussPos = xROI(miniGaussPosInd)
    
    smlROIind = find(xInt > max(xInt - pX(end))/2);
    smlROI = xROI(smlROIind);
    
    pXdbl = fminsearch(@(p)FitDoubleGauss(p,mnGaussPos,xROI(smlROIind),xInt(smlROIind)'), ...
        [pXtmp(1:2) 50 -10 3 min(xInt)],options) % position is passed to function separately so it's forced
    % to be at the correct place.
    pXdblComb = [pXdbl(1:4) mnGaussPos pXdbl(5:6)];
    
else
    pY = zeros(1,4);
    pX = pY;
end

figure(1)
subplot(3,3,[1 2 4 5])
oneMat = ones(size(OD));
zeroMat = nan(size(OD));
zeroMat(yROI,xROI) = 1;
ODwindow = OD.*zeroMat;
imagesc(ODwindow);
line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
axis equal
xlims = xlim;
ylims = ylim;

subplot(3,3,[3 6])
plot(yInt,yROI,gauss(pY,yC),yC);
set(gca,'YDir','reverse');
ylim(ylims);

subplot(3,3,[7 8])
plot(xROI,xInt,xC,gauss(pX,xC),xROI,DoubleGauss(pXdblComb,xROI));
xlim(xlims);

h = subplot(3,3,9,'Visible', 'off');
cla(h);
text(0, .5, sprintf('%s\n%s\n', ...
        'Mini-Gaussian Ampl.: ', ...
        num2str(pXdbl(4))), ...
        'VerticalAlignment', 'middle', ...
        'FontName', 'Courier New', ...
        'FontWeight', 'bold', ...
        'FontSize', 20);